library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.4     ✓ purrr   0.3.4
✓ tibble  3.1.2     ✓ dplyr   1.0.7
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.1
── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
library(SingleCellExperiment)
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats

Attaching package: ‘matrixStats’

The following object is masked from ‘package:dplyr’:

    count


Attaching package: ‘MatrixGenerics’

The following objects are masked from ‘package:matrixStats’:

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs,
    colLogSumExps, colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys,
    rowAvgsPerColSet, rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads,
    rowMaxs, rowMeans2, rowMedians, rowMins, rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs,
    rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians, rowWeightedSds, rowWeightedVars

Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply, parCapply, parLapply, parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect,
    is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
    table, tapply, union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:dplyr’:

    first, rename

The following object is masked from ‘package:tidyr’:

    expand

The following objects are masked from ‘package:base’:

    expand.grid, I, unname

Loading required package: IRanges

Attaching package: ‘IRanges’

The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice

The following object is masked from ‘package:purrr’:

    reduce

Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: ‘Biobase’

The following object is masked from ‘package:MatrixGenerics’:

    rowMedians

The following objects are masked from ‘package:matrixStats’:

    anyMissing, rowMedians
# Make sure that you use the forked version
# devtools::install_github("const-ae/sctransform@offset_with_flexible_theta")
library(sctransform)

Prepare Data

# Work around for some bug in zellkonverter
# https://github.com/theislab/zellkonverter/issues/45
setAs("dgRMatrix", to = "dgCMatrix", function(from){
  as(as(from, "CsparseMatrix"), "dgCMatrix")
})
if(! file.exists("../data/nih3t3.h5ad")){
  download.file("https://data.caltech.edu/tindfiles/serve/a448e98e-89cd-47b3-a134-803bbde29781/", "../data/nih3t3.h5ad")
}
if(! file.exists("../data/hek293t.h5ad")){
  download.file("https://data.caltech.edu/tindfiles/serve/b2758046-9247-43ab-b8f0-68882b4f39a3/", "../data/hek293t.h5ad")
}

Helper function, so that the content of the sctransform results stay in order

transform_sctrans_object <- function(sctrans_obj, Y){
  fit  <- list(overdispersions = rep(NA, nrow(Y)),
                  PearsonResiduals = array(NA, dim(Y)),
                  Mu <- array(NA, dim(Y)),
                  Beta = matrix(NA, nrow = nrow(Y), ncol = 2),
                  Beta_unreg = matrix(NA, nrow = nrow(Y), ncol = 2),
                  overdispersions_unreg = rep(NA, nrow(Y)),
                  gmean = rep(NA, nrow(Y)))
  sf <- colSums2(Y)
  
  indices <- match(rownames(sctrans_obj$model_pars_fit), rownames(Y))
  indices_unreg <- match(rownames(sctrans_obj$model_pars), rownames(Y))
  fit$overdispersions[indices] <- 1/sctrans_obj$model_pars_fit[,"theta"]
  fit$PearsonResiduals[indices, ] <- sctrans_obj$y
  fit$Beta[indices, ] <- sctrans_obj$model_pars_fit[, c("(Intercept)", "log_umi")]
  fit$Beta_unreg[indices_unreg, ] <- sctrans_obj$model_pars[, c("(Intercept)", "log_umi")]
  fit$Mu <- exp(fit$Beta %*% t(cbind(1, log10(sf))))
  fit$Mu_unreg <- exp(fit$Beta_unreg %*% t(cbind(1, log10(sf))))
  fit$overdispersions_unreg[indices_unreg] <- 1 / sctrans_obj$model_pars[, "theta"]
  fit$gmean[indices] <- sctrans_obj$gene_attr$gmean
  
  fit
}

Load Data

load_data_fnc <- list(
  NIH3T3 = function() zellkonverter::readH5AD("../data/nih3t3.h5ad"),
  HEK239T = function() zellkonverter::readH5AD("../data/hek293t.h5ad"),
  PBMC4k = function() TENxPBMCData::TENxPBMCData("pbmc4k"),
  Mesmer = function() scRNAseq::MessmerESCData(), 
  Zeisel = function() scRNAseq::ZeiselBrainData(),
  Zhen = function() DuoClustering2018::sce_filteredExpr10_Zhengmix4eq()
)

Run sctransform 3 times on each dataset


full_data <- bind_rows(lapply(names(load_data_fnc), function(name){
  print(paste0("Loading ", name))
  
  se <- load_data_fnc[[name]]()
  
  rownames(se) <- paste0("Gene_", 1:nrow(se))
  colnames(se) <- paste0("Cell_", 1:ncol(se))
  # Subsample data to 3000x1000 elements
  se_red <- se[sample(seq_len(nrow(se)), min(3000, nrow(se))), sample(seq_len(ncol(se)), min(1000, ncol(se)))]
  
  Y <- as.matrix(assay(se_red))
  
  set.seed(1)
  sctrans_intermediate <- vst(Y, method = "poisson")
  fit_pois <- transform_sctrans_object(sctrans_intermediate, Y)
  # To make sure that I am actually running the correct sctransform version
  # I check if the indicator value is set
  stopifnot(isTRUE(sctrans_intermediate$is_fork))
  
  set.seed(1)
  sctrans_intermediate2 <- vst(Y, method = "offset", theta_given = NULL)
  fit_off <- transform_sctrans_object(sctrans_intermediate2, Y)
  
  set.seed(1)
  sctrans_intermediate3 <- vst(Y, method = "offset", theta_given = 100)
  fit_off_fixed <- transform_sctrans_object(sctrans_intermediate3, Y)
  
  theta_given <- rep(100, nrow(Y))
  names(theta_given) <- rownames(Y)
  set.seed(1)
  sctrans_intermediate4 <- vst(Y, method = "nb_theta_given", theta_given = theta_given)
  fit_pois_fixed <- transform_sctrans_object(sctrans_intermediate4, Y)

  
  tibble(Dataset = name,
         Poisson = c(fit_pois$PearsonResiduals),
         PoissonFixed = c(fit_pois_fixed$PearsonResiduals),
         Offset = c(fit_off$PearsonResiduals),
         OffsetFixed = c(fit_off_fixed$PearsonResiduals),
         ExpressionLevel = c(fit_off$Mu),
         y = c(Y))
}))
[1] "Loading NIH3T3"
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 852 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 852 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
snapshotDate(): 2021-05-18

  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Found 41 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 852 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 13.36397 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 852 by 1000
Model formula is y ~ log_umi
Found 28 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 852 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 3.772232 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 852 by 1000
Model formula is y ~ log_umi
Second step: Get residuals using fitted parameters for 852 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 0.6714771 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 852 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 852 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Second step: Get residuals using fitted parameters for 852 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 4.452082 secs
[1] "Loading HEK239T"
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 927 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 927 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Found 38 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 927 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 4.839784 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 927 by 1000
Model formula is y ~ log_umi
Found 15 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 927 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 3.911496 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 927 by 1000
Model formula is y ~ log_umi
Second step: Get residuals using fitted parameters for 927 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 0.665138 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 927 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 927 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Found 1 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 927 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 4.321753 secs
[1] "Loading PBMC4k"
snapshotDate(): 2021-05-18
see ?TENxPBMCData and browseVignettes('TENxPBMCData') for documentation
loading from cache
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1127 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 1127 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===========================================================                                                                                                                     |  33%
  |                                                                                                                                                                                      
  |=====================================================================================================================                                                           |  67%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Found 80 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1127 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===========================================================                                                                                                                     |  33%
  |                                                                                                                                                                                      
  |=====================================================================================================================                                                           |  67%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 5.376468 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1127 by 1000
Model formula is y ~ log_umi
Found 49 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1127 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===========================================================                                                                                                                     |  33%
  |                                                                                                                                                                                      
  |=====================================================================================================================                                                           |  67%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 5.140258 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1127 by 1000
Model formula is y ~ log_umi
Second step: Get residuals using fitted parameters for 1127 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===========================================================                                                                                                                     |  33%
  |                                                                                                                                                                                      
  |=====================================================================================================================                                                           |  67%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 0.7097969 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1127 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 1127 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===========================================================                                                                                                                     |  33%
  |                                                                                                                                                                                      
  |=====================================================================================================================                                                           |  67%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Second step: Get residuals using fitted parameters for 1127 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===========================================================                                                                                                                     |  33%
  |                                                                                                                                                                                      
  |=====================================================================================================================                                                           |  67%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 5.405725 secs
[1] "Loading Mesmer"
snapshotDate(): 2021-05-18
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
snapshotDate(): 2021-05-18
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
snapshotDate(): 2021-05-18
loading from cache
require("ensembldb")
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1652 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 1652 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07
Found 27 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1652 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 5.554952 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1652 by 1000
Model formula is y ~ log_umi
Found 2 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1652 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 4.278342 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1652 by 1000
Model formula is y ~ log_umi
Second step: Get residuals using fitted parameters for 1652 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 0.7898989 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1652 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 1652 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Second step: Get residuals using fitted parameters for 1652 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 8.076696 secs
[1] "Loading Zeisel"
snapshotDate(): 2021-05-18
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
snapshotDate(): 2021-05-18
see ?scRNAseq and browseVignettes('scRNAseq') for documentation
loading from cache
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 2391 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
There are 1 estimated thetas smaller than 1e-07 - will be set to 1e-07
Found 4 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 2391 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===================================                                                                                                                                             |  20%
  |                                                                                                                                                                                      
  |======================================================================                                                                                                          |  40%
  |                                                                                                                                                                                      
  |==========================================================================================================                                                                      |  60%
  |                                                                                                                                                                                      
  |=============================================================================================================================================                                   |  80%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 7.448519 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 2391 by 1000
Model formula is y ~ log_umi
Found 3 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 2391 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===================================                                                                                                                                             |  20%
  |                                                                                                                                                                                      
  |======================================================================                                                                                                          |  40%
  |                                                                                                                                                                                      
  |==========================================================================================================                                                                      |  60%
  |                                                                                                                                                                                      
  |=============================================================================================================================================                                   |  80%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 7.491939 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 2391 by 1000
Model formula is y ~ log_umi
Found 2 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 2391 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===================================                                                                                                                                             |  20%
  |                                                                                                                                                                                      
  |======================================================================                                                                                                          |  40%
  |                                                                                                                                                                                      
  |==========================================================================================================                                                                      |  60%
  |                                                                                                                                                                                      
  |=============================================================================================================================================                                   |  80%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 0.9169271 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 2391 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 2000 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Found 1 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 2391 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |===================================                                                                                                                                             |  20%
  |                                                                                                                                                                                      
  |======================================================================                                                                                                          |  40%
  |                                                                                                                                                                                      
  |==========================================================================================================                                                                      |  60%
  |                                                                                                                                                                                      
  |=============================================================================================================================================                                   |  80%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 9.207349 secs
[1] "Loading Zhen"
see ?DuoClustering2018 and browseVignettes('DuoClustering2018') for documentation
loading from cache
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1556 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 1556 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Found 15 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1556 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 7.262152 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1556 by 1000
Model formula is y ~ log_umi
Found 18 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1556 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 7.017199 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1556 by 1000
Model formula is y ~ log_umi
Found 6 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1556 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 0.7863991 secs
Calculating cell attributes from input UMI matrix: log_umi
Variance stabilizing transformation of count matrix of size 1556 by 1000
Model formula is y ~ log_umi
Get Negative Binomial regression parameters per gene
Using 1556 genes, 1000 cells

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Found 2 outliers - those will be ignored in fitting/regularization step

Second step: Get residuals using fitted parameters for 1556 genes

  |                                                                                                                                                                                      
  |                                                                                                                                                                                |   0%
  |                                                                                                                                                                                      
  |============================================                                                                                                                                    |  25%
  |                                                                                                                                                                                      
  |========================================================================================                                                                                        |  50%
  |                                                                                                                                                                                      
  |====================================================================================================================================                                            |  75%
  |                                                                                                                                                                                      
  |================================================================================================================================================================================| 100%
Calculating gene attributes
Wall clock passed: Time difference of 7.148958 secs
p1 <- full_data %>%
  drop_na() %>%
  # sample_n(5000) %>%
  ggplot(aes(x= Poisson, y = Offset - Poisson)) +
    geom_hline(yintercept = 0, color = "#E9E9E9") +
    scattermore::geom_scattermore(aes(color = y), pointsize = 2.5, pixels = c(512, 512), alpha = 0.1) +
    scale_color_viridis_c(direction = 1, option = "magma", begin = 0.2, end = 0.8, 
                          trans = "log1p", breaks = c(0, 10, 100, 1000),
                          labels = c(0, 10, 100, ">1000"), limits = c(0, 1001), 
                          oob = scales::oob_squish_any, 
                          minor_breaks = c(0:9, seq(10, 90, by = 10), seq(100, 1000, by = 100))) +
    scale_y_continuous(limits = c(-30, 30)) +
    facet_wrap(~ Dataset, nrow = 1) +
    cowplot::theme_cowplot() +
    labs(title = expression(Effect~of~fixing~beta[s]==1),
         color = "Y", 
         x = "Residuals from the default sctransform model",
         y = expression(Delta[1]))


p1

p2 <- full_data %>%
  drop_na() %>%
  # sample_n(5000) %>%
  ggplot(aes(x= Poisson, y = PoissonFixed - Poisson)) +
    geom_hline(yintercept = 0, color = "#E9E9E9") +
    scattermore::geom_scattermore(aes(color = y), pointsize = 2.5, pixels = c(512, 512), alpha = 0.1) +
    scale_color_viridis_c(direction = 1, option = "magma", begin = 0.2, end = 0.8, 
                          trans = "log1p", breaks = c(0, 10, 100, 1000),
                          labels = c(0, 10, 100, ">1000"), limits = c(0, 1001), 
                          oob = scales::oob_squish_any, 
                          minor_breaks = c(0:9, seq(10, 90, by = 10), seq(100, 1000, by = 100))) +
    scale_y_continuous(limits = c(-30, 30)) +
    facet_wrap(~ Dataset, nrow = 1) +
    cowplot::theme_cowplot() +
    labs(title = expression(Effect~of~fixing~alpha==0.01),
         color = "Y", 
         x = "Residuals from the default sctransform model",
         y = expression(Delta[2]))

p2

p3 <- full_data %>%
  drop_na() %>%
  # sample_n(5000) %>%
  ggplot(aes(x= Poisson, y = OffsetFixed - Poisson)) +
    geom_hline(yintercept = 0, color = "#E9E9E9") +
    scattermore::geom_scattermore(aes(color = y), pointsize = 2.5, pixels = c(512, 512), alpha = 0.1) +
    scale_color_viridis_c(direction = 1, option = "magma", begin = 0.2, end = 0.8, 
                          trans = "log1p", breaks = c(0, 10, 100, 1000),
                          labels = c(0, 10, 100, ">1000"), limits = c(0, 1001), 
                          oob = scales::oob_squish_any, 
                          minor_breaks = c(0:9, seq(10, 90, by = 10), seq(100, 1000, by = 100))) +
    scale_y_continuous(limits = c(-30, 30)) +
    facet_wrap(~ Dataset, nrow = 1) +
    cowplot::theme_cowplot() +
    labs(title = expression(Effect~of~fixing~beta[s]==1~and~alpha==0.01),
         color = "Y", 
         x = "Residuals from the default sctransform model",
         y = expression(Delta[3]))

p3

library(patchwork)
res <- (p1 / p2 / p3) +
  plot_annotation(tag_levels = "A") + plot_layout(guides = "collect")

res

cowplot::save_plot("../output/sctransform_beta_s_and_alpha_importance.pdf", res, 
                   ncol = 7, nrow = 3, base_asp = 0.6, base_height = 3, device = cairo_pdf)

Session Info

sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] patchwork_1.1.1             DuoClustering2018_1.10.0    ensembldb_2.16.0            AnnotationFilter_1.16.0     GenomicFeatures_1.44.0      AnnotationDbi_1.54.1       
 [7] scRNAseq_2.6.1              TENxPBMCData_1.10.0         HDF5Array_1.20.0            rhdf5_2.36.0                DelayedArray_0.18.0         Matrix_1.3-4               
[13] sctransform_0.3.2.9002      SingleCellExperiment_1.14.1 SummarizedExperiment_1.22.0 Biobase_2.52.0              GenomicRanges_1.44.0        GenomeInfoDb_1.28.0        
[19] IRanges_2.26.0              S4Vectors_0.30.0            BiocGenerics_0.38.0         MatrixGenerics_1.4.0        matrixStats_0.59.0          forcats_0.5.1              
[25] stringr_1.4.0               dplyr_1.0.7                 purrr_0.3.4                 readr_1.4.0                 tidyr_1.1.3                 tibble_3.1.2               
[31] ggplot2_3.3.4               tidyverse_1.3.1            

loaded via a namespace (and not attached):
  [1] readxl_1.3.1                  backports_1.2.1               AnnotationHub_3.0.0           BiocFileCache_2.0.0           plyr_1.8.6                    lazyeval_0.2.2               
  [7] BiocParallel_1.26.0           listenv_0.8.0                 scattermore_0.7               digest_0.6.27                 htmltools_0.5.1.1             viridis_0.6.1                
 [13] fansi_0.5.0                   magrittr_2.0.1                memoise_2.0.0                 globals_0.14.0                Biostrings_2.60.1             modelr_0.1.8                 
 [19] prettyunits_1.1.1             colorspace_2.0-1              blob_1.2.1                    rvest_1.0.0                   rappdirs_0.3.3                haven_2.4.1                  
 [25] xfun_0.24                     crayon_1.4.1                  RCurl_1.98-1.3                jsonlite_1.7.2                glue_1.4.2                    gtable_0.3.0                 
 [31] zlibbioc_1.38.0               XVector_0.32.0                Rhdf5lib_1.14.1               future.apply_1.7.0            scales_1.1.1                  DBI_1.1.1                    
 [37] ggthemes_4.2.4                Rcpp_1.0.6                    viridisLite_0.4.0             xtable_1.8-4                  progress_1.2.2                reticulate_1.20              
 [43] mclust_5.4.7                  bit_4.0.4                     httr_1.4.2                    dir.expiry_1.0.0              ellipsis_0.3.2                farver_2.1.0                 
 [49] pkgconfig_2.0.3               XML_3.99-0.6                  dbplyr_2.1.1                  utf8_1.2.1                    labeling_0.4.2                tidyselect_1.1.1             
 [55] rlang_0.4.11                  reshape2_1.4.4                later_1.2.0                   munsell_0.5.0                 BiocVersion_3.13.1            cellranger_1.1.0             
 [61] tools_4.1.0                   cachem_1.0.5                  cli_2.5.0                     generics_0.1.0                RSQLite_2.2.7                 ExperimentHub_2.0.0          
 [67] broom_0.7.7                   fastmap_1.1.0                 yaml_2.2.1                    knitr_1.33                    bit64_4.0.5                   fs_1.5.0                     
 [73] KEGGREST_1.32.0               future_1.21.0                 mime_0.10                     xml2_1.3.2                    biomaRt_2.48.1                compiler_4.1.0               
 [79] rstudioapi_0.13               filelock_1.0.2                curl_4.3.1                    png_0.1-7                     interactiveDisplayBase_1.30.0 reprex_2.0.0                 
 [85] stringi_1.6.2                 basilisk.utils_1.4.0          lattice_0.20-44               ProtGenerics_1.24.0           vctrs_0.3.8                   pillar_1.6.1                 
 [91] lifecycle_1.0.0               rhdf5filters_1.4.0            BiocManager_1.30.16           cowplot_1.1.1                 bitops_1.0-7                  httpuv_1.6.1                 
 [97] rtracklayer_1.52.0            R6_2.5.0                      BiocIO_1.2.0                  promises_1.2.0.1              gridExtra_2.3                 parallelly_1.26.0            
[103] codetools_0.2-18              zellkonverter_1.2.0           MASS_7.3-54                   assertthat_0.2.1              rjson_0.2.20                  withr_2.4.2                  
[109] GenomicAlignments_1.28.0      Rsamtools_2.8.0               GenomeInfoDbData_1.2.6        hms_1.1.0                     grid_4.1.0                    basilisk_1.4.0               
[115] shiny_1.6.0                   lubridate_1.7.10              restfulr_0.0.13              
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShTaW5nbGVDZWxsRXhwZXJpbWVudCkKIyBNYWtlIHN1cmUgdGhhdCB5b3UgdXNlIHRoZSBmb3JrZWQgdmVyc2lvbgojIGRldnRvb2xzOjppbnN0YWxsX2dpdGh1YigiY29uc3QtYWUvc2N0cmFuc2Zvcm1Ab2Zmc2V0X3dpdGhfZmxleGlibGVfdGhldGEiKQpsaWJyYXJ5KHNjdHJhbnNmb3JtKQpgYGAKCiMgUHJlcGFyZSBEYXRhCgpgYGB7cn0KIyBXb3JrIGFyb3VuZCBmb3Igc29tZSBidWcgaW4gemVsbGtvbnZlcnRlcgojIGh0dHBzOi8vZ2l0aHViLmNvbS90aGVpc2xhYi96ZWxsa29udmVydGVyL2lzc3Vlcy80NQpzZXRBcygiZGdSTWF0cml4IiwgdG8gPSAiZGdDTWF0cml4IiwgZnVuY3Rpb24oZnJvbSl7CiAgYXMoYXMoZnJvbSwgIkNzcGFyc2VNYXRyaXgiKSwgImRnQ01hdHJpeCIpCn0pCmBgYAoKYGBge3J9CmlmKCEgZmlsZS5leGlzdHMoIi4uL2RhdGEvbmloM3QzLmg1YWQiKSl7CiAgZG93bmxvYWQuZmlsZSgiaHR0cHM6Ly9kYXRhLmNhbHRlY2guZWR1L3RpbmRmaWxlcy9zZXJ2ZS9hNDQ4ZTk4ZS04OWNkLTQ3YjMtYTEzNC04MDNiYmRlMjk3ODEvIiwgIi4uL2RhdGEvbmloM3QzLmg1YWQiKQp9CmlmKCEgZmlsZS5leGlzdHMoIi4uL2RhdGEvaGVrMjkzdC5oNWFkIikpewogIGRvd25sb2FkLmZpbGUoImh0dHBzOi8vZGF0YS5jYWx0ZWNoLmVkdS90aW5kZmlsZXMvc2VydmUvYjI3NTgwNDYtOTI0Ny00M2FiLWI4ZjAtNjg4ODJiNGYzOWEzLyIsICIuLi9kYXRhL2hlazI5M3QuaDVhZCIpCn0KYGBgCgoKCkhlbHBlciBmdW5jdGlvbiwgc28gdGhhdCB0aGUgY29udGVudCBvZiB0aGUgc2N0cmFuc2Zvcm0gcmVzdWx0cyBzdGF5IGluIG9yZGVyCmBgYHtyfQp0cmFuc2Zvcm1fc2N0cmFuc19vYmplY3QgPC0gZnVuY3Rpb24oc2N0cmFuc19vYmosIFkpewogIGZpdCAgPC0gbGlzdChvdmVyZGlzcGVyc2lvbnMgPSByZXAoTkEsIG5yb3coWSkpLAogICAgICAgICAgICAgICAgICBQZWFyc29uUmVzaWR1YWxzID0gYXJyYXkoTkEsIGRpbShZKSksCiAgICAgICAgICAgICAgICAgIE11IDwtIGFycmF5KE5BLCBkaW0oWSkpLAogICAgICAgICAgICAgICAgICBCZXRhID0gbWF0cml4KE5BLCBucm93ID0gbnJvdyhZKSwgbmNvbCA9IDIpLAogICAgICAgICAgICAgICAgICBCZXRhX3VucmVnID0gbWF0cml4KE5BLCBucm93ID0gbnJvdyhZKSwgbmNvbCA9IDIpLAogICAgICAgICAgICAgICAgICBvdmVyZGlzcGVyc2lvbnNfdW5yZWcgPSByZXAoTkEsIG5yb3coWSkpLAogICAgICAgICAgICAgICAgICBnbWVhbiA9IHJlcChOQSwgbnJvdyhZKSkpCiAgc2YgPC0gY29sU3VtczIoWSkKICAKICBpbmRpY2VzIDwtIG1hdGNoKHJvd25hbWVzKHNjdHJhbnNfb2JqJG1vZGVsX3BhcnNfZml0KSwgcm93bmFtZXMoWSkpCiAgaW5kaWNlc191bnJlZyA8LSBtYXRjaChyb3duYW1lcyhzY3RyYW5zX29iaiRtb2RlbF9wYXJzKSwgcm93bmFtZXMoWSkpCiAgZml0JG92ZXJkaXNwZXJzaW9uc1tpbmRpY2VzXSA8LSAxL3NjdHJhbnNfb2JqJG1vZGVsX3BhcnNfZml0WywidGhldGEiXQogIGZpdCRQZWFyc29uUmVzaWR1YWxzW2luZGljZXMsIF0gPC0gc2N0cmFuc19vYmokeQogIGZpdCRCZXRhW2luZGljZXMsIF0gPC0gc2N0cmFuc19vYmokbW9kZWxfcGFyc19maXRbLCBjKCIoSW50ZXJjZXB0KSIsICJsb2dfdW1pIildCiAgZml0JEJldGFfdW5yZWdbaW5kaWNlc191bnJlZywgXSA8LSBzY3RyYW5zX29iaiRtb2RlbF9wYXJzWywgYygiKEludGVyY2VwdCkiLCAibG9nX3VtaSIpXQogIGZpdCRNdSA8LSBleHAoZml0JEJldGEgJSolIHQoY2JpbmQoMSwgbG9nMTAoc2YpKSkpCiAgZml0JE11X3VucmVnIDwtIGV4cChmaXQkQmV0YV91bnJlZyAlKiUgdChjYmluZCgxLCBsb2cxMChzZikpKSkKICBmaXQkb3ZlcmRpc3BlcnNpb25zX3VucmVnW2luZGljZXNfdW5yZWddIDwtIDEgLyBzY3RyYW5zX29iaiRtb2RlbF9wYXJzWywgInRoZXRhIl0KICBmaXQkZ21lYW5baW5kaWNlc10gPC0gc2N0cmFuc19vYmokZ2VuZV9hdHRyJGdtZWFuCiAgCiAgZml0Cn0KYGBgCgoKTG9hZCBEYXRhCgpgYGB7cn0KbG9hZF9kYXRhX2ZuYyA8LSBsaXN0KAogIE5JSDNUMyA9IGZ1bmN0aW9uKCkgemVsbGtvbnZlcnRlcjo6cmVhZEg1QUQoIi4uL2RhdGEvbmloM3QzLmg1YWQiKSwKICBIRUsyMzlUID0gZnVuY3Rpb24oKSB6ZWxsa29udmVydGVyOjpyZWFkSDVBRCgiLi4vZGF0YS9oZWsyOTN0Lmg1YWQiKSwKICBQQk1DNGsgPSBmdW5jdGlvbigpIFRFTnhQQk1DRGF0YTo6VEVOeFBCTUNEYXRhKCJwYm1jNGsiKSwKICBNZXNtZXIgPSBmdW5jdGlvbigpIHNjUk5Bc2VxOjpNZXNzbWVyRVNDRGF0YSgpLCAKICBaZWlzZWwgPSBmdW5jdGlvbigpIHNjUk5Bc2VxOjpaZWlzZWxCcmFpbkRhdGEoKSwKICBaaGVuID0gZnVuY3Rpb24oKSBEdW9DbHVzdGVyaW5nMjAxODo6c2NlX2ZpbHRlcmVkRXhwcjEwX1poZW5nbWl4NGVxKCkKKQpgYGAKClJ1biBzY3RyYW5zZm9ybSAzIHRpbWVzIG9uIGVhY2ggZGF0YXNldAoKLSBTdGFuZGFyZCBQb2lzc29uIG1vZGVsCi0gT2Zmc2V0IHdpdGggZmxleGlibGUgdGhldGEgZXN0aW1hdGlvbgotIE9mZnNldCB3aXRoIGZpeGVkIHRoZXRhCgpgYGB7ciwgd2FybmluZyA9IEZBTFNFfQoKZnVsbF9kYXRhIDwtIGJpbmRfcm93cyhsYXBwbHkobmFtZXMobG9hZF9kYXRhX2ZuYyksIGZ1bmN0aW9uKG5hbWUpewogIHByaW50KHBhc3RlMCgiTG9hZGluZyAiLCBuYW1lKSkKICAKICBzZSA8LSBsb2FkX2RhdGFfZm5jW1tuYW1lXV0oKQogIAogIHJvd25hbWVzKHNlKSA8LSBwYXN0ZTAoIkdlbmVfIiwgMTpucm93KHNlKSkKICBjb2xuYW1lcyhzZSkgPC0gcGFzdGUwKCJDZWxsXyIsIDE6bmNvbChzZSkpCiAgIyBTdWJzYW1wbGUgZGF0YSB0byAzMDAweDEwMDAgZWxlbWVudHMKICBzZV9yZWQgPC0gc2Vbc2FtcGxlKHNlcV9sZW4obnJvdyhzZSkpLCBtaW4oMzAwMCwgbnJvdyhzZSkpKSwgc2FtcGxlKHNlcV9sZW4obmNvbChzZSkpLCBtaW4oMTAwMCwgbmNvbChzZSkpKV0KICAKICBZIDwtIGFzLm1hdHJpeChhc3NheShzZV9yZWQpKQogIAogIHNldC5zZWVkKDEpCiAgc2N0cmFuc19pbnRlcm1lZGlhdGUgPC0gdnN0KFksIG1ldGhvZCA9ICJwb2lzc29uIikKICBmaXRfcG9pcyA8LSB0cmFuc2Zvcm1fc2N0cmFuc19vYmplY3Qoc2N0cmFuc19pbnRlcm1lZGlhdGUsIFkpCiAgIyBUbyBtYWtlIHN1cmUgdGhhdCBJIGFtIGFjdHVhbGx5IHJ1bm5pbmcgdGhlIGNvcnJlY3Qgc2N0cmFuc2Zvcm0gdmVyc2lvbgogICMgSSBjaGVjayBpZiB0aGUgaW5kaWNhdG9yIHZhbHVlIGlzIHNldAogIHN0b3BpZm5vdChpc1RSVUUoc2N0cmFuc19pbnRlcm1lZGlhdGUkaXNfZm9yaykpCiAgCiAgc2V0LnNlZWQoMSkKICBzY3RyYW5zX2ludGVybWVkaWF0ZTIgPC0gdnN0KFksIG1ldGhvZCA9ICJvZmZzZXQiLCB0aGV0YV9naXZlbiA9IE5VTEwpCiAgZml0X29mZiA8LSB0cmFuc2Zvcm1fc2N0cmFuc19vYmplY3Qoc2N0cmFuc19pbnRlcm1lZGlhdGUyLCBZKQogIAogIHNldC5zZWVkKDEpCiAgc2N0cmFuc19pbnRlcm1lZGlhdGUzIDwtIHZzdChZLCBtZXRob2QgPSAib2Zmc2V0IiwgdGhldGFfZ2l2ZW4gPSAxMDApCiAgZml0X29mZl9maXhlZCA8LSB0cmFuc2Zvcm1fc2N0cmFuc19vYmplY3Qoc2N0cmFuc19pbnRlcm1lZGlhdGUzLCBZKQogIAogIHRoZXRhX2dpdmVuIDwtIHJlcCgxMDAsIG5yb3coWSkpCiAgbmFtZXModGhldGFfZ2l2ZW4pIDwtIHJvd25hbWVzKFkpCiAgc2V0LnNlZWQoMSkKICBzY3RyYW5zX2ludGVybWVkaWF0ZTQgPC0gdnN0KFksIG1ldGhvZCA9ICJuYl90aGV0YV9naXZlbiIsIHRoZXRhX2dpdmVuID0gdGhldGFfZ2l2ZW4pCiAgZml0X3BvaXNfZml4ZWQgPC0gdHJhbnNmb3JtX3NjdHJhbnNfb2JqZWN0KHNjdHJhbnNfaW50ZXJtZWRpYXRlNCwgWSkKCiAgCiAgdGliYmxlKERhdGFzZXQgPSBuYW1lLAogICAgICAgICBQb2lzc29uID0gYyhmaXRfcG9pcyRQZWFyc29uUmVzaWR1YWxzKSwKICAgICAgICAgUG9pc3NvbkZpeGVkID0gYyhmaXRfcG9pc19maXhlZCRQZWFyc29uUmVzaWR1YWxzKSwKICAgICAgICAgT2Zmc2V0ID0gYyhmaXRfb2ZmJFBlYXJzb25SZXNpZHVhbHMpLAogICAgICAgICBPZmZzZXRGaXhlZCA9IGMoZml0X29mZl9maXhlZCRQZWFyc29uUmVzaWR1YWxzKSwKICAgICAgICAgRXhwcmVzc2lvbkxldmVsID0gYyhmaXRfb2ZmJE11KSwKICAgICAgICAgeSA9IGMoWSkpCn0pKQoKCmBgYAoKCgoKCgpgYGB7ciwgZmlnLndpZHRoPTE1LCBmaWcuaGVpZ2h0PTN9CnAxIDwtIGZ1bGxfZGF0YSAlPiUKICBkcm9wX25hKCkgJT4lCiAgIyBzYW1wbGVfbig1MDAwKSAlPiUKICBnZ3Bsb3QoYWVzKHg9IFBvaXNzb24sIHkgPSBPZmZzZXQgLSBQb2lzc29uKSkgKwogICAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMCwgY29sb3IgPSAiI0U5RTlFOSIpICsKICAgIHNjYXR0ZXJtb3JlOjpnZW9tX3NjYXR0ZXJtb3JlKGFlcyhjb2xvciA9IHkpLCBwb2ludHNpemUgPSAyLjUsIHBpeGVscyA9IGMoNTEyLCA1MTIpLCBhbHBoYSA9IDAuMSkgKwogICAgc2NhbGVfY29sb3JfdmlyaWRpc19jKGRpcmVjdGlvbiA9IDEsIG9wdGlvbiA9ICJtYWdtYSIsIGJlZ2luID0gMC4yLCBlbmQgPSAwLjgsIAogICAgICAgICAgICAgICAgICAgICAgICAgIHRyYW5zID0gImxvZzFwIiwgYnJlYWtzID0gYygwLCAxMCwgMTAwLCAxMDAwKSwKICAgICAgICAgICAgICAgICAgICAgICAgICBsYWJlbHMgPSBjKDAsIDEwLCAxMDAsICI+MTAwMCIpLCBsaW1pdHMgPSBjKDAsIDEwMDEpLCAKICAgICAgICAgICAgICAgICAgICAgICAgICBvb2IgPSBzY2FsZXM6Om9vYl9zcXVpc2hfYW55LCAKICAgICAgICAgICAgICAgICAgICAgICAgICBtaW5vcl9icmVha3MgPSBjKDA6OSwgc2VxKDEwLCA5MCwgYnkgPSAxMCksIHNlcSgxMDAsIDEwMDAsIGJ5ID0gMTAwKSkpICsKICAgIHNjYWxlX3lfY29udGludW91cyhsaW1pdHMgPSBjKC0zMCwgMzApKSArCiAgICBmYWNldF93cmFwKH4gRGF0YXNldCwgbnJvdyA9IDEpICsKICAgIGNvd3Bsb3Q6OnRoZW1lX2Nvd3Bsb3QoKSArCiAgICBsYWJzKHRpdGxlID0gZXhwcmVzc2lvbihFZmZlY3R+b2Z+Zml4aW5nfmJldGFbc109PTEpLAogICAgICAgICBjb2xvciA9ICJZIiwgCiAgICAgICAgIHggPSAiUmVzaWR1YWxzIGZyb20gdGhlIGRlZmF1bHQgc2N0cmFuc2Zvcm0gbW9kZWwiLAogICAgICAgICB5ID0gZXhwcmVzc2lvbihEZWx0YVsxXSkpCgoKcDEKYGBgCgpgYGB7ciwgZmlnLndpZHRoPTE1LCBmaWcuaGVpZ2h0PTN9CnAyIDwtIGZ1bGxfZGF0YSAlPiUKICBkcm9wX25hKCkgJT4lCiAgIyBzYW1wbGVfbig1MDAwKSAlPiUKICBnZ3Bsb3QoYWVzKHg9IFBvaXNzb24sIHkgPSBQb2lzc29uRml4ZWQgLSBQb2lzc29uKSkgKwogICAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gMCwgY29sb3IgPSAiI0U5RTlFOSIpICsKICAgIHNjYXR0ZXJtb3JlOjpnZW9tX3NjYXR0ZXJtb3JlKGFlcyhjb2xvciA9IHkpLCBwb2ludHNpemUgPSAyLjUsIHBpeGVscyA9IGMoNTEyLCA1MTIpLCBhbHBoYSA9IDAuMSkgKwogICAgc2NhbGVfY29sb3JfdmlyaWRpc19jKGRpcmVjdGlvbiA9IDEsIG9wdGlvbiA9ICJtYWdtYSIsIGJlZ2luID0gMC4yLCBlbmQgPSAwLjgsIAogICAgICAgICAgICAgICAgICAgICAgICAgIHRyYW5zID0gImxvZzFwIiwgYnJlYWtzID0gYygwLCAxMCwgMTAwLCAxMDAwKSwKICAgICAgICAgICAgICAgICAgICAgICAgICBsYWJlbHMgPSBjKDAsIDEwLCAxMDAsICI+MTAwMCIpLCBsaW1pdHMgPSBjKDAsIDEwMDEpLCAKICAgICAgICAgICAgICAgICAgICAgICAgICBvb2IgPSBzY2FsZXM6Om9vYl9zcXVpc2hfYW55LCAKICAgICAgICAgICAgICAgICAgICAgICAgICBtaW5vcl9icmVha3MgPSBjKDA6OSwgc2VxKDEwLCA5MCwgYnkgPSAxMCksIHNlcSgxMDAsIDEwMDAsIGJ5ID0gMTAwKSkpICsKICAgIHNjYWxlX3lfY29udGludW91cyhsaW1pdHMgPSBjKC0zMCwgMzApKSArCiAgICBmYWNldF93cmFwKH4gRGF0YXNldCwgbnJvdyA9IDEpICsKICAgIGNvd3Bsb3Q6OnRoZW1lX2Nvd3Bsb3QoKSArCiAgICBsYWJzKHRpdGxlID0gZXhwcmVzc2lvbihFZmZlY3R+b2Z+Zml4aW5nfmFscGhhPT0wLjAxKSwKICAgICAgICAgY29sb3IgPSAiWSIsIAogICAgICAgICB4ID0gIlJlc2lkdWFscyBmcm9tIHRoZSBkZWZhdWx0IHNjdHJhbnNmb3JtIG1vZGVsIiwKICAgICAgICAgeSA9IGV4cHJlc3Npb24oRGVsdGFbMl0pKQoKcDIKYGBgCgoKCmBgYHtyLCBmaWcud2lkdGg9MTUsIGZpZy5oZWlnaHQ9M30KcDMgPC0gZnVsbF9kYXRhICU+JQogIGRyb3BfbmEoKSAlPiUKICAjIHNhbXBsZV9uKDUwMDApICU+JQogIGdncGxvdChhZXMoeD0gUG9pc3NvbiwgeSA9IE9mZnNldEZpeGVkIC0gUG9pc3NvbikpICsKICAgIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDAsIGNvbG9yID0gIiNFOUU5RTkiKSArCiAgICBzY2F0dGVybW9yZTo6Z2VvbV9zY2F0dGVybW9yZShhZXMoY29sb3IgPSB5KSwgcG9pbnRzaXplID0gMi41LCBwaXhlbHMgPSBjKDUxMiwgNTEyKSwgYWxwaGEgPSAwLjEpICsKICAgIHNjYWxlX2NvbG9yX3ZpcmlkaXNfYyhkaXJlY3Rpb24gPSAxLCBvcHRpb24gPSAibWFnbWEiLCBiZWdpbiA9IDAuMiwgZW5kID0gMC44LCAKICAgICAgICAgICAgICAgICAgICAgICAgICB0cmFucyA9ICJsb2cxcCIsIGJyZWFrcyA9IGMoMCwgMTAsIDEwMCwgMTAwMCksCiAgICAgICAgICAgICAgICAgICAgICAgICAgbGFiZWxzID0gYygwLCAxMCwgMTAwLCAiPjEwMDAiKSwgbGltaXRzID0gYygwLCAxMDAxKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgb29iID0gc2NhbGVzOjpvb2Jfc3F1aXNoX2FueSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgbWlub3JfYnJlYWtzID0gYygwOjksIHNlcSgxMCwgOTAsIGJ5ID0gMTApLCBzZXEoMTAwLCAxMDAwLCBieSA9IDEwMCkpKSArCiAgICBzY2FsZV95X2NvbnRpbnVvdXMobGltaXRzID0gYygtMzAsIDMwKSkgKwogICAgZmFjZXRfd3JhcCh+IERhdGFzZXQsIG5yb3cgPSAxKSArCiAgICBjb3dwbG90Ojp0aGVtZV9jb3dwbG90KCkgKwogICAgbGFicyh0aXRsZSA9IGV4cHJlc3Npb24oRWZmZWN0fm9mfmZpeGluZ35iZXRhW3NdPT0xfmFuZH5hbHBoYT09MC4wMSksCiAgICAgICAgIGNvbG9yID0gIlkiLCAKICAgICAgICAgeCA9ICJSZXNpZHVhbHMgZnJvbSB0aGUgZGVmYXVsdCBzY3RyYW5zZm9ybSBtb2RlbCIsCiAgICAgICAgIHkgPSBleHByZXNzaW9uKERlbHRhWzNdKSkKCnAzCmBgYAoKCgpgYGB7ciwgZmlnLndpZHRoPTIwLCBmaWcuaGVpZ2h0PTEwfQpsaWJyYXJ5KHBhdGNod29yaykKcmVzIDwtIChwMSAvIHAyIC8gcDMpICsKICBwbG90X2Fubm90YXRpb24odGFnX2xldmVscyA9ICJBIikgKyBwbG90X2xheW91dChndWlkZXMgPSAiY29sbGVjdCIpCgpyZXMKYGBgCgoKYGBge3J9CmNvd3Bsb3Q6OnNhdmVfcGxvdCgiLi4vb3V0cHV0L3NjdHJhbnNmb3JtX2JldGFfc19hbmRfYWxwaGFfaW1wb3J0YW5jZS5wZGYiLCByZXMsIAogICAgICAgICAgICAgICAgICAgbmNvbCA9IDcsIG5yb3cgPSAzLCBiYXNlX2FzcCA9IDAuNiwgYmFzZV9oZWlnaHQgPSAzLCBkZXZpY2UgPSBjYWlyb19wZGYpCmBgYAoKCgoKCiMgU2Vzc2lvbiBJbmZvCgpgYGB7cn0Kc2Vzc2lvbkluZm8oKQpgYGAKCgog